
cd('.')
addpath('matlab functions')

D0 = readtable('NC Data/SCC_4to5.csv');

for mdl = {'math_NOmatching','math_matching','reading_NOmatching','reading_matching'}

    vamdl = load(sprintf('output/main/vam_%s.mat',mdl{1}));
    modspecs = vamdl.modspecs;
    tchcoef = modspecs.tchcoef;
    D = D0(~isnan(D0.(modspecs.outcome)),:);
    vadata = genvadata(D,modspecs,vamdl.controls);

    vahat = cell([numel(vadata) 1]);
    for j = 1:numel(vadata)
        vahat{j} = zeros([numel(vadata(j).Y) numel(tchcoef) 2]);
        [~,sid_rw] = ismember(vadata(j).sid,D.sid);
        yr = D.year(sid_rw);
        for yr_val = unique(yr)'
            this_year = (yr==yr_val);
            last_year = (yr==(yr_val-1));
            next_year = (yr==(yr_val+1));
            for m = 1:2
                if m==1
                    I = ~(this_year | last_year);
                elseif m==2
                    I = ~(this_year | next_year) ;
                end
                vad = vadata(j);
                vad.Y = vad.Y(I);
                vad.Z = vad.Z(I,:);
                vad.ZF = vad.ZF(I,:);
                vad.sid = vad.sid(I);
                EB = vaposterior(vamdl.VAdistEstimates,vad);
                EB_dev = EB - vamdl.VAmomentEstimates.mean;
                vahat{j}(this_year,:,m) = ones(nnz(this_year),1)*EB_dev';        
            end
        end
    end

    vahat = cat(1,vahat{:});
    vahat = reshape(vahat,size(vahat,1),[]);
    OUT = array2table([cat(1,vadata.sid) vahat], ...
        'VariableNames',[{'sid'} ...
                         strcat('coef_2yr_l_',tchcoef) ...
                         strcat('coef_2yr_f_',tchcoef)]);
    OUT = innerjoin(D(:,[{'sid','year','districtid','schoolid','teacherid','schgrade'} ...
        modspecs.outcome tchcoef]), OUT, 'Keys', 'sid');
    OUT = renamevars(OUT,modspecs.outcome, 'score');
    OUT.tv_2yr_l = sum(OUT{:,strcat('coef_2yr_l_',tchcoef)}.*OUT{:,tchcoef},2);
    OUT.tv_2yr_f = sum(OUT{:,strcat('coef_2yr_f_',tchcoef)}.*OUT{:,tchcoef},2);
    OUT = groupsummary(OUT,{'districtid','schoolid','year','schgrade'},'mean',{'score','tv_2yr_l','tv_2yr_f'});
    OUT = renamevars(OUT,{'GroupCount','mean_tv_2yr_f','mean_tv_2yr_l','mean_score'},{'n','tv_2yr_f','tv_2yr_l','score'});

    writetable(OUT,sprintf('cfr test/tv_%s.csv',mdl{1}))

end

